%每次随机翻动flip_num个格点，再用Metropolis算法对体系状态变化进行选择
function [mat, k, m] = trans_state(mat, n, T, b, flip_num)
    J = 1;
    mat0 = mat;
    %E0 = ising_energy(mat, n, T, b);
    i=0;
    num = round(flip_num);
    %将循环次数上限设置为2 * n * n /num，循环结束后若还未找到合适状态，则依旧使用原状态
    %这样能防止在能量局域最小构型下的长时间循环
    while i < (2 * n * n / num)
        i=i+1;
        selection = ceil(n * n * rand(1,num));
        mid = mat(selection);
        mat(selection) = -mid;
        [k,m] = ind2sub(size(mat),selection);
        dE = sum(- 2 * b * (-mid)/ T - 2 * J / T * (-mid) .* near_sum(mat, n, k, m));
        if(dE < 0)
            return
        else
            rnum = rand;
            if(rnum < exp(-dE))
                return
            else
                mat = mat0;
            end
        end
    end
    mat = mat0;
    k = zeros(1,num) ;
    m = zeros(1,num) ;
    return
end

function sumi = near_sum(mat, n, k, m)
    len = length(k);
    sumi = zeros(1,len);
    for x = 1:len
        if(k(x) == 1)
            if(m(x)==1)
                sumi(x) = mat(1, 2) + mat(2, 1) + mat(1, n) + mat(n, 1);
            elseif(m(x)==n)
                sumi(x) = mat(1, n-1) + mat(2, n) + mat(n, n) + mat(1, 1);
            else
                sumi(x) = mat(1, m(x)-1) + mat(1, m(x)+1) + mat(2, m(x)) + mat(n, m(x));
            end
        elseif(k(x) == n)
            if(m(x)==1)
                sumi(x) = mat(n, 2) + mat(n-1, 1) + mat(1, 1) + mat(n, n);
            elseif(m(x)==n)
                sumi(x) = mat(n, n-1) + mat(n-1, n) + mat(1, n) + mat(n, 1);
            else
                sumi(x) = mat(n, m(x)-1) + mat(n, m(x)+1) + mat(n-1, m(x)) + mat(1, m(x));
            end
        else
            if(m(x)==1)
                sumi(x) = mat(k(x)-1, 1) + mat(k(x)+1, 1) + mat(k(x), 2) + mat(k(x), n);
            elseif(m(x)==n)
                sumi(x) = mat(k(x)-1, n) + mat(k(x)+1, n) + mat(k(x), n-1) + mat(k(x), 1);
            else
                sumi(x) = mat(k(x)-1, m(x)) + mat(k(x)+1, m(x)) + mat(k(x), m(x)-1) + mat(k(x), m(x)+1);
            end
        end
    end
    return
end
